Surface diffusion coefficients by thermodynamic integration: Cu on Cu(100) 
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The rate of diffusion of a Cu adatom on the Cu(100) surface is calculated using thermodynamic 
integration within the transition state theory. The results are found to be in excellent agreement 
with the essentially exact values from molecular-dynamics simulations. The activation energy and 
related entropy are shown to be effectively independent of temperature, thus establishing the validity 
of the Arrhenius law over a wide range of temperatures. Our study demonstrates the equivalence of 
diffusion rates calculated using thermodynamic integration within the transition state theory and 
direct molecular-dynamics simulations. 
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Precise knowledge of diffusion processes is essential to 
understanding non-equilibrium phenomena such as nu- 
cleation and growth [Q. On surfaces, for instance, the 
rates at which particles diffuse determine the equilib- 
rium shape of islands and, on macroscopic timescales, 
the morphology of films. Yet, very little is known of 
the fundamentals of diffusion. Diffusion constants, for 
one, are notably difficult to measure and accurate data 
is available only for the simplest mechanisms on a small 
number of simple surfaces Because diffusion is an 
activated (Arrhenius) process (at low enough tempera- 
tures Q), small errors in the energy barriers translate 
into large uncertainties in the diffusion coefficients, and 
thus surface structure. In addition, in order to deter- 
mine the pre-exponential factor, several measurements 
are needed in a range of temperatures over which the 
Arrhenius behaviour is expected to hold, which is not 
always feasible: in practice, the value of the prefactor is 
often prescribed. This is a dangerous state of affairs since 
diffusion obeys the Meyer-Neldel compensation law jl| - 
for a family of related processes, the prefactor increases 
exponentially with the activation barrier ^ . 

On the theory side, the situation is just as difficult. It 
is necessary, in order to describe diffusion accurately, to 
have a proper model for the interatomic potentials. Semi- 
empirical models, such as the embedded-atom method 
(EAM) Q, while simple and sometimes remarkably ac- 
curate, lack the transferability and predictive power of 
first-principles methods. The latter, however, are sub- 
ject to size and other limitations, and uncertainties are 
difficult to estimate. For instance, even for such a simple 
case as diffusion by jumps of Cu adatoms on the Cu(100) 
surface, experiment and ab initio calculations disagree 
0; the origin of the discrepancy remains unclear. 



Because of various limitations, the technique used for 
computing diffusion rates also is important. The sim- 
plest option consists in simulating diffusion explicitely 
using molecular dynamics (MD). One advantage of this 
method is that a priori knowledge of the diffusion mech- 
anism is not required. Such calculations are however too 
demanding for ab initio methods. Also, the simulations 
have to be carried out at relatively high temperatures 
where diffusion is "active" on MD timescales; at high 
temperature, however, diffusion often proceeds by the 
combination of several mechanisms, making it difficult 
to extract individual contributions. Finally, because of 
possible anharmonic contributions, the calculated Arrhe- 
nius law may not extrapolate to low temperatures. 

Another option consists in computing directly the acti- 
vation barrier and the prefactor using the transition-state 
theory (TST) and various approximations p-pT|. Here, 
however, the reaction path must be known; while this 
might be a limitation for bulk diffusion, it is usually not 
a serious problem for surfaces where diffusion is relatively 
well characterized. In the context of TST, and given a 
model for the interatomic potentials, free-energy calcula- 
tions, in particular thermodynamic integration (TI) , offer 
the most accurate route to the study of diffusion pro- 
cesses. In this approach, the diffusion path is followed 
step by step, and the free energy calculated using finite- 
temperature MD. The procedure works best at low tem- 
perature; at high temperature, indeed, diffusion events 
are more frequent and the atoms must be constrained 
to their equilibrium positions (see below). In this case, 
it might be more advantageous to use the explicit MD 
approach. 

Because the two methods are so different, and cover 
different temperature ranges, and because diffusion is an 



1 



important, difficult, and yet unresolved problem in most 
cases, it is of utmost interest to ascertain that they lead 
to equivalent results. This question has been addressed 
previously using Monte Carlo simulations with restricted 
dynamics on Lennard- Jones metals OJl^], but the re- 
sults were not conclusive: the energy barriers were found 
to differ by as much as 35% and the prefactors by a 
factor of ~1.8. Here we reexamine the problem in the 
case of Cu diffusion on Cu(100), for which detailed MD 
simulations with EAM potentials have recently been re- 
ported 0; EAM provides a rather accurate description 
of the energetics of Cu. The TI is performed in full us- 
ing MD, solving directly the TST equations. We find the 
explicit MD and the TST/TI calculations to be in very 
close agreement for both the prefactor and the energy 
barrier. The free-energy barrier, in addition, is found to 
depend linearly on temperature, confirming the validity 
of the Arrhenius law over a wide range of temperatures. 
Our results establish unambiguously the equivalence of 
the two methods, thus providing a useful framework for 
the calculation of diffusion constants. 

In the TST, the rate of reaction from one equilibrium 
site to another, via a saddle point, is given by Jl5| ] 

k = K-k TST , k TST = ve- AW / kBT , (1) 

where k is the transmission coefficient (or "recrossing 
rate") and &tst is the TST rate constant. AW is the 
activation free energy; the prefactor v, the frequency at 
which the reaction is attempted, is given by 

(2) 
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The integral in Eq. (Q) runs between two transition sites 
a distance L apart, say from Xb — L to Xb, via the equilib- 
rium site at x m . W(x) is the "potential of mean force": 
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where < /(A) > is the mean force that must be applied 
in order to constrain the particle at position A along the 
reaction path; evidently < / > is zero if x = x m or 
x = Xb- W can be obtained numerically by calculating 
the mean force at several points along the diffusion path 
using constrained MD Q. 

TABLE I. Comparison between TI and MD results for 
the jump (J) and exchange (X) diffusion activation barriers 
AE (in eV) and rate prefactors To (in THz); also given are 
the entropy AS (in fcs) and the static energy barrier, AE(0). 
Estimated errors are given in parenthesis. 
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FIG. 1. (a) Activation free energy vs temperature for 
jumps (squares, dashed line) and exchanges (circles, full line); 
the lines are linear fits to the finite-temperature points, (b) 
Attempt-to-diffuse frequencies vs temperature; the lines are 
the predictions of the simple model discussed in the text, (c) 
Transmission coefficients vs temperature. 

The TI calculations were carried out using MD and 
EAM potentials. As in Ref. |Q, the surface was modeled 
by a slab consisting of eight layers, each containing 64 
atoms, with the bottom two fixed in their equilibrium 
lattice positions; periodic boundary conditions were ap- 
plied in the two lateral directions. We investigate here 
the four temperatures 100, 300, 500, and 800 K; this will 
permit a comparison with our earlier MD calculations, 
which covered the range 650-850 K . Most calculations 
were done in the NVT ensemble, using a Nose thermo- 
stat to control the temperature JTq] ; however, we have 
also done some calculations in the NVE ensemble to as- 
sess the effect of the thermostat. At each point along 
the reaction path, the system was first equilibrated for 
48 ps, then statistics accumulated for a further f20 ps. 
At the highest temperatures, the atoms lying close to 
that undergoing diffusion were attached to their equilib- 
rium positions with harmonic springs. Several values of 
the spring constant were examined and the mean force 
obtained by extrapolating to zero (!(]] . 

The transmission coefficient k is given by 



K =< Q[x(+t) - X b ] - Q[x(-t) - X b ] >i»r vib 



(4) 



where T v ;b is a time characteristic of atomic vibrations 
and O is the Heaviside step function, k was obtained by 
averaging over f 00 different initial configurations, taken 
at 1.2 ps intervals from a MD run with the adatom con- 
strained at the saddle point. Each of these was run for 
f .2 ps both backward and forward in time p7| . 

We plot in Fig. |l](a) the activation free energies AW as 
a function of temperature for both mechanisms possible 
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FIG. 2. Diffusion rates vs inverse temperature; the lines 
are fits to an Arrhenius law; estimated errors are smaller than 
the size of the symbols. Inset: Comparison between MD and 
TST/TI results in the temperature range relevant to the MD 
data. 

on this surface, viz., jump and exchange; the static (0 
K) values are also indicated. In both cases, AW is very 
well represented by a linear function of temperature, i.e., 
AW = AE - TAS, where AE and AS are both, ef- 
fectively, temperature- independent. The values of AE 
and AS are listed in Table | along with the correspond- 
ing values for AE from the direct MD simulations. For 
both diffusion processes, the TI and MD barriers are in 
excellent agreement with the static barriers. 

We note from Fig. 0(a) that, in spite of the fact that it 
has a larger activation barrier, exchange diffusion is more 
favourable than jump diffusion above 800 K or so. This is 
a manifestation of the Meyer- Neldel rule [0,|| ; the prefac- 
tor for exchanges is much larger (20 times — cf. Table |) 
than that for jumps, thus compensating for the smaller 
activation term. Compensation is so efficient that the 
process with a larger barrier becomes dominant at suf- 
ficiently high temperature. Thus, even at low tempera- 
ture, where anharmonic effects are small, multi-phononic 
contributions to the entropy cannot be ignored. 

The attempt-to-diffuse frequencies for the two pro- 
cesses are displayed in Fig. 0(b). They depend only very 
weakly on temperature. The deviations at 800 K might 
reflect some error in the TI calculations at high tempera- 
tures. However, even taking such possible errors into ac- 
count, the observed temperature dependence is insignifi- 
cant relative to the exponential activation term and thus, 
for all practical purposes, the prefactor for diffusion via 
a given mechanism can be taken as constant. 

The (slow) variation of v with temperature can be un- 
derstood in terms of the following simple anharmonic 
model: Taking W(x) to be of the form ax 2 /2 — fix 3 /3, 
with W{x m ) = and AW = W(x b ) - W(x m ) (so that 



a = 6AW/xl and [3 = a>AW/x\) one easily finds from 
Eq. (0) (neglecting the anharmonic term in the evalua- 
tion of the integral) that, in the low temperature limit, 
v(T) = u Q y/AE- TAS/VAE with u = v(0). We see 
from Fig. 0(b) that the TI data is very well fitted by 
this model up to about 500 K; as expected, this ap- 
proximation is no longer valid at higher temperature. 
Note that the slight temperature dependence will be neg- 
ligible, again, on an Arrhenius plot. The differences 
in the ^-values for the two processes arise, to a large 
extent, from "geometrical" differences. For the above 
model we also have 2Ttv = mxi, where m is 

the mass of the diffusing entity — mc u for jump and 
tocu/2 for exchange (motion of a dimer with respect to 
its center of mass); taking x b ,x = 1.6 A for an exchange 
(roughly a/2, with a = 3.61 A the lattice parameter) 
and Xb j — 1.3 A for a jump [a/(2\/2~)], we find, indeed, 
AWj/xIj « AW x /xl x . " 

The transmission coefficient is the probability that a 
diffusion event actually takes place once the saddle point 
is reached. For both mechanisms, k depends relatively 
little on temperature, as can be seen in Fig. 0(c). For 
jumps, the transmission coefficient is about 0.9, and thus 
has little effect on the diffusion prefactor. For exchanges, 
k is close to 0.6, and the effect is slightly more important. 

The diffusion rate is the product of transition rate, 
transmission coefficient and number of equivalent reac- 
tion paths. (The diffusion constant is obtained from the 
diffusion rate by multiplying by a geometrical factor.) 
For both jumps and exchanges, there are four equivalent 
paths and we thus have [cf. Eq. (ED)]: 

r = Anue-^ w ' kBT = T e~ AE / kBT ,T = 4^e AS / feB (5) 

The TI results for T are presented in Fig. 0. The data are 
extremely well fitted by an Arrhenius law at all tempera- 
tures, even as large as 800 K. The resulting values of AE 
are nearly identical to those determined earlier by fitting 
to the free energies. The slight temperature dependence 
of the attempt-to-diffuse frequencies and the transmis- 
sion coefficients has, as anticipated, no visible effect on 
the Arrhenius barriers. The values of the prefactors r , 
which we return to below, are listed in Table |[ 

There has been some concerns that the thermostat in 
NVT simulations might lead to sizable errors in free- 
energy calculations (see, e.g., Ref. fig]). In order to test 
this, we have carried out some TI calculations for both 
jump and exchange at 500 K, using both NVT and NVE 
algorithms. Differences were found to be insignificant 
- at most 0.007 eV on free energies and 0.01 THz on 
diffusion rates — well within numerical uncertainties. 

In the inset of Fig. 0, finally, we compare closely the 
TI results with the MD simulations. The former covers 
the range 0-800 K, while the latter is for 650-850 K. 
The TI and MD calculations are found to be in complete 
agreement for both diffusion mechanisms over the whole 
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temperature range. This establishes without ambiguity 
that the two different computational schemes comple- 
ment one another exactly. In addition, our calculations 
demonstrate that the range of validity of the Arrhenius 
law can extend over a much wider range of temperatures 
than is normally assumed. 

Free-energy calculations of the barriers for jump diffu- 
sion on Cu and Ag (100) surfaces based on the harmonic 
approximation to TST have been reported recently [^0) . 
The calculations were carried out using the same model 
potentials (EAM) as in the present study; yet, for Cu 
jump on Cu(100), a prefactor 10 times smaller than that 
found here was obtained. Numerical error cannot be to- 
tally excluded as the cause for this discrepancy, but the 
consistency between our TI and MD results strongly sug- 
gests that this is not the case. Rather, it is more likely 
a problem with methodology: the harmonic and quasi- 
harmonic approximations neglect the multiphononic con- 
tributions which affect deeply the thermodynamic func- 
tions, especially prefactors, giving rise, as we have seen 
earlier, to such effects as the Meyer- Neldel law Q. 

It has been claimed by many authors (see for example 
Refs. [0, H] and 10) that the entropy AS and the 
energy barrier AE depend on temperature. Our results 
provide no evidence for this. The separation of the differ- 
ent terms in Eq. ([!]) is somewhat arbitrary and largely 
a matter of definition. The simplest expression for T, 
viz. r = T cxp(-AE/k B T), where T (and thus AS) as 
well as AE are effectively independent of temperature, is 
able to account very precisely for both the TI and the 
MD data over the full range of temperatures considered. 
Indeed, the entropy term, after dividing by ksT, merely 
renormalizes the prefactor [cf. Eq. (0)]. 

We have reported a detailed comparison of the rates 
for jump and exchange self-diffusion on Cu(100) as ob- 
tained from full thermodynamic integration and direct 
molecular-dynamics simulations. We find the two meth- 
ods to be in perfect agreement over a wide range of tem- 
peratures. Our results clearly demonstrate that a sim- 
ple representation of the diffusion rate in terms of a 
static energy barrier (which defines the activation term) 
and a temperature-independent entropy (which defines 
the prefactor), as they appear in the usual transition 
state theory, accounts fully for the dynamics of isolated 
adatoms. Furthermore, the present study clearly demon- 
strates the equivalence of the diffusion constants obtained 
within TST/TI and from direct MD simulations. 
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